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We compare three different methods to determine the lattice spacing in lattice QCD and give 
results from calculations on the MILC ensembles of configurations that include the effect of u, d 
and s sea quarks. It is useful, for ensemble to ensemble comparison, to express the results as giving 
a physical value for ri , a parameter from the heavy quark potential. Combining the three methods 
gives a value for ri in the continuum limit of 0.3133(23) (3) fm. Using the MILC values for ro/r\, 
this corresponds to a value for the ro parameter of 0.4661(38) fm. We also discuss how to use the 
rj s for determining the lattice spacing and tuning the s-quark mass accurately, by giving values for 
m„ s (0.6858(40) GeV) and /„ s (0.1815(10) GeV). 



I. INTRODUCTION 

Results from lattice QCD calculations are generally 
computed in units of the lattice spacing a used in the 
simulation. The lattice spacing must be computed sep- 
arately and divided out in order to convert these results 
into physical units (GeV, fm...), for comparison with 
experiment. Any error in the lattice spacing determina- 
tion feeds into most other quantities from lattice QCD, 
and, in many cases, it is among the dominant sources 
of errors. For example, in our determination of the de- 
cay constant of the D s meson p], 1% of the total error of 
1.3% comes from the 1.5% uncertainty in the value of the 
lattice spacing. Reducing the error on the lattice spacing 
is then very important for increasing the precision of the 
realistic lattice QCD calculations now possible j2]. 

Generally the value of the lattice spacing is determined 
by comparing values from the simulation, in lattice units, 
with values from experiment, in physical units. A lattice 
simulation, for example, might give a value for the pion 
decay constant in lattice units: a/^ at . Dividing by the ex- 
perimental value /° xp in GeV gives a value for the lattice 
spacing, a — (a/^. at )//° xp , in inverse GeV. This lattice 
spacing can then be used to convert other simulation re- 
sults from lattice units to physical units. 

Lattice spacings determined in this way are inherently 
ambiguous because lattice simulations are never exact. 
In particular the use of a nonzero lattice spacing causes 
lattice quantites, like /^. at , to deviate from their physical 
values, in this case /^ xp . Such errors differ from quantity 
to quantity, and therefore so will values for the lattice 
spacing that are computed from these quantities. Such 
differences, however, vanish in the continuum limit, a — > 
0, and so do not affect lattice predictions that have been 
extrapolated to a = 0. 
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In principle, any dimensionful quantity can be used to 
determine the lattice spacing, but some quantities are 
more useful than others. Ideally one wants quantities 
that are easily computed, free of other types of simu- 
lation error, largely independent of lattice parameters 
other than the lattice spacing, and well measured in ex- 
periments. Use of the pion decay constant, for example, 
is not ideal. This decay constant is quite sensitive to the 
u and d quark masses, which are generally too large in 
current simulations; accurate values for the decay con- 
stant can be obtained only after chiral extrapolations of 
the simulation data to the physical quark masses. This 
greatly complicates the use of the decay constant to set 
the lattice spacing. 

One physical quantity that is very easy to calculate 
in lattice simulations is the r\ parameter derived from 
the potential V{r) between two infinite-mass quarks sep- 
arated by distance r. Parameter r\ is defined implicitly 
by the equation rlF(ri) — C where F(r) = dV/dr and 
C = 1 |3J. (Taking C = 1.65 gives the original such 
standard parameter, r [4..) This quantity is easily cal- 
culated, in lattice units (that is, ri/a), to better than 1%. 
Unlike the pion decay constant, it is only weakly depen- 
dent upon the quark masses. It would be an ideal choice 
for setting the lattice spacing except for the fact that 
there is no experimental value for the physical ri — this 
must be estimated instead from other lattice calculations. 

In this paper we examine three other quantities that 
can be used to determine the lattice spacing: 1) the 
radial excitation energy in the T system (toy< — mx); 
2) the mass difference between the D s meson and one 
half the T] c mass; and 3) the decay constant of the ficti- 
tious rj s particle, which can be related accurately to fx 
and / w . The valence-quark masses are easily tuned in 
each case and each quantity is relatively insensitive to 
sea-quark masses. Consequently each of these quantities 
can be used to generate lattice spacings on an ensemble- 
by-ensemble basis. 

None of these quantities can be computed as accurately 
as r\j a in simulations, but we can combine simulation re- 
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suits for them with values for r\/ato obtain very accurate 
estimates for the physical value of r\ . Given r\ , the differ- 
ent values of r\/a can be used to obtain accurate lattice 
spacings for each of the simulations we discuss here and 
any other simulations where r\ja has been computed. 

Of our three quantities, the r\ s decay constant gives 
the most accurate results. The r\ s is a fictitious meson, 
however, and so its "experimental" properties must be 
related to those of real mesons using simulations. The 
i] s is particularly closely related to the ir and K mesons. 
As we will show, its mass and decay constant can be 
accurately related to those of the tt and K through a 
chiral analysis of simulation data for a variety of quark 
masses and lattice spacings. Such an analysis also gives 
an independent, fourth estimate of r%. 

We describe in section 2 the three primary methods we 
have used to obtain lattice spacings for a wide variety of 
simulations. Each can be used to generate an estimate 
for the physical value of r 1; given values of r^ja. In 
section 3 we combine the three analyses to generate a 
single, combined estimate for r\. This can then be used 
to covert the r\ja values into a determination of a on 
each ensemble. We also demonstrate how to determine 
the lattice spacing from the rj s without using r-y . The two 
methods are compared and shown to agree in the a — > 
limit. In Section 4 we give a value for ro derived from our 
value of r\ for comparison to others using that parameter. 
In section 5 summarize our results. Finally, we discuss 
the chiral analysis of decay constants and masses for the 
7r, K and their relation with those of the r] s meson in 
Appendices A, B and C. 



II. LATTICE CALCULATION 

In Table [I] we list the parameter sets for the different 
MILC ensembles of gluon configurations that we have 
used here, although not all ensembles were used in every 
lattice spacing determination. 

Values for the static-quark potential parameter ri/a, 
in lattice units, were determined by the MILC collabo- 
ration [5] . They calculated the heavy quark potential by 
fitting Wilson loops of fixed spatial size as a function of 
lattice time. On the finest two sets of ensembles smeared 
time links were used to reduce statistical noise and a two- 
state exponential fit in time reduced the contamination 
from excited potentials. The heavy quark potential ob- 
tained was then fit as a function of spatial separation over 
the range between 0.2 fm and 0.7 fm to a Cornell poten- 
tial with the addition of corrections for lattice artifacts. 
The point at which the condition for r\ held was then 
determined from this fit. The errors given are statistical 
errors only, since discretisation effects are taken care of 
in our continuum extrapolations. 

In what follows we will combine these values for r\/a 
with estimates of the lattice spacing a determined using 
three different physical quantities to obtain estimates for 
the physical value of r\ (that is, at zero lattice spacing 



and with correct sea-quark masses) . 



A. m T / — mr 

The calculation of the spectrum of mesons formed as 
bound states of bottom quarks and antiquarks has been 
an important test for lattice QCD. There are many radial 
and orbital excitations below threshold for strong decay 
and so many gold-plated states, well-characterised ex- 
perimentally The radial and orbital excitation energies 
are almost identical for charmonium and bottomonium 
when spin-averaged [B] and so rather insensitive to the 
heavy quark mass. Heavy-quark vacuum polarization ef- 
fects are tiny and so can be safely neglected. This makes 
these systems very suitable for the determination of the 
lattice spacing [7| and was one of the key calculations 
demonstrating the importance of including the effect of 
u, d and s sea quarks [2]. 

Here we improve on the calculations in [7 which used 
results from MILC super-coarse, coarse and fine en- 
sembles and compared ensembles with and without sea 
quarks. We study only ensembles including sea quarks 
but include also very coarse and superfine ensembles for 
a wider range of lattice spacing values. 

We calculate 6-quark propagators on the MILC gluon 
field configurations using lattice NonRelativistic QCD 
(NRQCD) which has been developed over many years 
to handle well the physics of heavy quark systems on the 
lattice [8]. It makes a virtue of the nonrelativistic nature 
of bottomonium bound states (i> 2 ss 0.1 for the T) by 
discarding the rest mass energy in favour of accurately 
handling typical momentum and energy scales inside the 
bound states. NRQCD can be matched to full QCD or- 
der by order in v 2 and a s . We work through 0(v^) in 
the nonrelativistic expansion and apply discretisation im- 
provements through 0(a 2 ) to v\ terms and to chromo- 
magnetic and chromoelectric field-dependent terms at vf 
(so that terms which induce fine structure are completely 
improved to C(a 4 )). An analysis of remaining systematic 
errors is given in |7]. 

The NonRelativistic Hamiltonian that we use is given 
by 0: 

aH = aHo + a5H; 
A( 2 ) 

aH -- 



2aM h 



(A< 2 >) 2 ig 



8(oJkf t ) 3 8(<zM 6 ) 2 



V-E - E- V 



-c 3 



8(aM b )- 



V x E - E x V 



9 ^ a 2 A (4) 
-C4- — — a ■ B + C5; 



-C 6 : 



2aM b 2AaM b 
a (A( 2 )) 2 



(1) 



'Wn{aM b ) 2 ' 

This is implemented in calculating b quark propagators 
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TABLE I: Ensembles (sets) of MILC configurations with 
gauge coupling f3, size L a xT and sea mass parameters mf aq 
and mf q used for this analysis. The sea ASQTAD quark 
masses (I — u/d) are given in the MILC convention where 
uo is the plaquette tadpole parameter. The lattice spacing 
values in units of ri after 'smoothing' are given in the third 
column [5]. Sets 1 and 2 are 'very coarse'; sets 3, 4 and 5, 
'coarse'; sets 6 and 7 'fine'; set 8 'superfine' and set 9 'ultra- 
fine'. 



Set 


P 


ri/a 


asq 

auom l 


au mf q 


L/a 


T/a 


1 


6.572 


2.152(5) 


0.0097 


0.0484 


16 


48 


2 


6.586 


2.138(4) 


0.0194 


0.0484 


16 


48 


3 


6.76 


2.647(3) 


0.005 


0.05 


24 


64 


4 


6.76 


2.618(3) 


0.01 


0.05 


20 


64 


5 


6.79 


2.644(3) 


0.02 


0.05 


20 


64 


6 


7.09 


3.699(3) 


0.0062 


0.031 


28 


96 


7 


7.11 


3.712(4) 


0.0124 


0.031 


28 


96 


8 


7.46 


5.296(7) 


0.0036 


0.018 


48 


144 


9 


7.81 


7.115(20) 


0.0028 


0.014 


64 


192 



by evolving the 6 quark Green's function on a single pass 
through the lattice using: 

(i-^r(i-^)G(sr,*) (2) 



with starting condition: 

G(x,0) = 4>(x)l. 



(3) 



Here V is the symmetric lattice derivative and V is the 

A( 2 > is the 



improved derivative, Vfc = — Vjj. 3 '/6. 
standard lattice discretisation of the second derivative 
J2 3 Vf 5 and AW is J2j V f } - aM b is the bare b quark 
mass in lattice units. 

4>{x) is a real spatial smearing function which multi- 
plies a unit matrix in color and spin space as the starting 
point for the quark propagator. The antiquark propa- 
gator for a given source is then the complex conjugate 
of the Green's function obtained from eq. [2j When the 
quark and antiquark propagators are combined (with ap- 
propriate Pauli matrices for different J PC [10]) into me- 
son correlators, <p improves the overlap with particular 
ground and excited states for a better signal. This will 
be discussed further below. 

n is a stability parameter which is chosen to tame (un- 
physical) high momentum modes of the b quark propaga- 
tor which might otherwise cause the meson correlators to 
grow exponentially with time rather than fall. We used 
n = 4 throughout this calculation instead of the value 
ri = 2 used earlier [7J, so that we could work on finer 
lattices and keep the same value of n for all ensembles. 
This means that discretisation errors are smoothly con- 
nected from one lattice spacing to another and the higher 
value of n has the advantage of reducing some systematic 
errors. 



TABLE II: Parameters used in our calculations of b quark 
propagators and bb correlators on various MILC ensembles, 
numbered as in Table [T] Mf,a is the bare b quark mass in 
lattice units, uol is the tadpole-improvement factor and the 
stability factor n is taken as 4 everywhere. n c / 9 gives the 
number of gluon field configurations used from the ensemble 
and n t is the number of time sources for b quark propagators 
per configuration. T is the time length in lattice units of the 
propagators, ao is the size parameter for the quark smearing 
function (f> es (x) given in eq. [4] 



Set 


aM b 


uol 




n t 


T 


a 


1 


3.4 


0.8218 


631 


24 


32 


0.83 


2 


3.4 


0.8225 


631 


24 


32 


0.83 


3 


2.8 


0.8362 


2083 


32 


32 


1.0 


4 


2.8 


0.8359 


595 


32 


32 


1.0 


6 


1.95 


0.8541 


557 


8 


48 


1.41 


8 


1.34 


0.8696 


698 


8 


48 


2.0 



We use 'tadpole-improvement' [TT] for all terms by di- 
viding the gluon fields by a factor uq when they are 
read in. uq is taken as uql, the value of the mean trace of 
the gluon field for that ensemble in lattice Landau gauge 
(where the trace is maximised). This removes, in a mean- 
field way, the disparity between lattice and continuum 
gluon fields induced by the fact that the lattice field is 
exponentially related to the continuum field. Single com- 
posite operators, such as A^ 4 \ are expanded out fully 
so that all cancellations between U and IF are correctly 
tadpole-improved. This is not done for the (l — aH /2n) n 
terms. Table [TT] gives our parameters for the ensembles 
used in the T analysis. 

In order to reduce statistical errors over our previous 
calculation we have investigated a number of improve- 
ments. The first was to look at different forms for the 
quark smearing <fi(x). The simplest is a S function but in 
addition we can take an arbitrary functional form for <p(x) 
provided that the gluon field configurations are gauge- 
fixed, at least on a time-slice. The MILC configurations 
that we use here are fixed to Coulomb gauge. When a 
b quark propagator from a S source and a b propagator 
from a <f> = f(x) source are combined a good overlap with 
a particular T state is expected when, in the language of 
a potential model, f(x) is a good approximation to the 
wavefunction of that state. The ground state T(1S) will 
dominate all 1 correlators eventually so that there is 
no advantage in including a smearing function that gives 
a good overlap with that state [12]. Instead, to obtain a 
good signal for the 25 — IS splitting, we concentrated on 
functions that had very small overlap with the IS state, 
and therefore had better information about radial exci- 
tations. A very good smearing for this was the function 
from [7J called (p es : 



4> es (r) = (2a - r) exp(-r/(2a )). 



(4) 



The size parameter, ao, was tuned on coarse lattices to 
reduce the overlap of the correlator (known as the 'ee' 
correlator, see below) with the ground state, as judged 
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by the small amplitude of the correlator at large times 
when the ground state dominates, ao was then scaled 
as appropriate to ensembles of different lattice spacing. 
Values are given in Table [nj By combining b and b prop- 
agators from 8 function sources and <p es sources we are 
able to make up 3 different meson smearing functions: 
I is from combining two 8 sources; e is from combining 
a <p es source with a 8 source and E is from combining 
two (j) es sources (so that the composite meson smearing 
function is then the convolution of <f> es with itself). I, e 
and E smearing functions can also be applied at the sink 
to make a 3 x 3 matrix of correlators, with notation 11, 
le, ee etc. 

We also used a random wall source for our b quark 
propagators, taking a set of U(l) random numbers, r, 
with unit norm at every point on a time slice, one set 
for each color of the b quark propagator. These were 
combined with the smearing functions <f> so that 

G(x,0) ClCl = ^2(f)(\x - y\)r(cx,y)l apin . (5) 
y 

When quark and antiquark propagators are combined to- 
gether the random noise cancels except where the initial 
spatial sites are the same and this effectively increases the 
number of meson correlators sampled. We find that the 
error on the ground state T energy is reduced by a factor 
of 3 on coarse lattices and 5 on fine lattices, when cor- 
rected to the same number of configurations. The excited 
state energy does not improve by the same factor, how- 
ever. Indeed we found rather little improvement in the 
error on excited state energies which mirrors our experi- 
ence with applying random wall sources to B mesons |13j . 
The inference is that random wall sources are much less 
effective in situations where the degradation of the sig- 
nal/noise is exponential. We calculate propagators from 
many different time sources (which we then average over) 
per configuration to improve statistical precision further. 
The details of numbers of configurations and time sources 
are collected in Table [ill 

As in [7] we use a Bayesian fitting method [14] to fit the 
3x3 matrix of hadron correlators to a multi-exponential 
form to extract the energies of states appearing in that 
correlator. This alows us to fit the entire correlator (i.e. 
for all time separations between source and sink) , so mak- 
ing use of all the information contained in it. It also 
means that the fit results we obtain, for example for 
the ground state, include the effect of the higher excited 
states that are present in the correlator, and are not bi- 
assed by an attempt to fit only one or two states in a 
particular time window. The fitting function is 

G meson (n sc , n sk ;t) = ^ a(n sc , k)a*(n sk , k)e~ Ekt . (6) 
fe=i 

where a(n sc / sk , k) are the (real) amplitudes for state k 
to appear in the smearings used at the source and sink 
of the correlator respectively. 



The Bayesian fitting method [T3] allows a large num- 
ber of exponentials to be used in the fit by constraining 
the way in which these exponentials can appear based on 
physical information. The simplest physical information 
is that the energies of states are ordered, and we imple- 
ment this in the fit by taking the energy fit parameters 
as the natural logarithms of the ground state energy and 
of the energy splittings between adjacent states. On top 
of this we apply priors to the splittings between adjacent 
states that constrain them to be of order 500 MeV with 
a width of a factor of two, i.e. between 250 MeV and 
1000 MeV. Amplitudes are typically constrained around 
zero with a width of 1.0 (our composite meson smearing 
functions are normalised so that the spatial sum of their 
square is 1). We apply a cut on the range of eigenvalues 
present in the correlation matrix of 10 -3 except for the 
high statistics calculation on the coarse 005/05 lattices 
where we use 10~ 4 . This reduces the number of degrees 
of freedom in the fit to between 120 and 170, with 208 in 
the coarse 005/05 fit. We obtain values for the T ground 
state energy and that of the first radial excitation, the T' 
as a function of the number of exponentials in the fits. 
We demand a good x 2 and that the fit for 3 adjacent 
exponentials should agree both on the fitted values for 
the energies of interest and on the errors. The ground 
state energy stabilises very quickly, but the first excited 
state is not generally stable until we reach 8 exponen- 
tials. Fit results on the different MILC ensembles are 
then tabulated in Table [ITT] from 10 exponential fits. 

Figure[l]shows results from our highest statistics calcu- 
lation on the coarse 005/05 lattices. Here we are able to 
obtain a good signal for even higher excited states than 
the 25. The plot shows the ratio of the 35 — 15 splitting 
and the 45 - 15 splitting to that of the 25 - 15. The 
35 — 15 splitting is obtained to 3% and in agreement 
with experiment. The 45— 15 splitting is not very accu- 
rate even with the statistics available here. The result is 
slightly higher than experiment, but the 45 state is not 
gold-plated, decaying to BB. This is not taken account 
of accurately in the lattice calculation and so we expect 
our result to be higher than experiment. In our lower 
statistics calculations we do not generally have a signifi- 
cant signal for the 45 and our 35—15 splitting has an 
error of between 5% and 10%. 

As discussed earlier, the excitation energies for bound 
states of heavy quarks are almost independent of the 
heavy quark mass, meaning that accurate tuning of this 
mass is not required for these splittings. Use of the ran- 
dom wall does, however, allow us to determine the meson 
energy as a function of meson momentum much more ac- 
curately than in previous calculations, and so the meson 
'kinetic' masses can be well determined. The meson mass 
in NRQCD must be determined from the meson disper- 
sion relation because the zero of energy has an offset. 
The mass is then given by the difference in energy be- 
tween mesons at zero momentum and momentum pa on 
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FIG. 1: Results for highly excited states from our fit to the 
3x3 matrix of T correlators from the coarse 005/05 (set 4) 
ensemble as a function of the number of exponentials included 
in the fit. The /dof is also shown - the fit had 208 degrees 
of freedom. The results are stable from 9 to 12 exponentials. 



the lattice by [7]: 



Ma = 



p 2 a 2 - (AEa) 2 
2AEa ' 



(7) 



AEa is calculated by taking the difference in energy of 
the ground state from a simultaneous fit to 11 T corre- 
lators made from a standard (zero momentum) random 
wall as described above, and a random wall patterned 
with an appropriate Fourier factor to give a momentum 
of (1,0,0) to both quark and antiquark, so that the T 
has momentum (2,0,0) (and its permutations). Values 
obtained for the kinetic mass on each of the MILC en- 
sembles are given in Table |III| They are tuned within 
10% of t he e xperimental result of 9.46 GeV [15] . 

Table III gives results for the lattice spacings a T ob- 
tained by dividing the simulation results for a(E 2 — Ei) 
by the experimental value of 0.5630 GeV for the split- 
ting. Statistical errors are at the level of 1%. System- 
atic errors arise from two sources, discretisation errors 
and missing higher-order relativistic correctiosn to the 
NRQCD action. The former can be removed by con- 
tinuum extrapolation as long as they are well-behaved. 
The leading discretisation corrections come from radia- 
tive corrections to existing terms in the action and can be 
calculated in perturbation theory. They have been shown 
to be small corrections in the region of aMb in which we 
work, and relatively independent of aM& [16] . Relativis- 
tic corrections survive the continuum limit and are the 
main source of systematic error for this method. They 



TABLE III: Results for the ground state energy, aEi, and 
radial excitation energy, a(E2 — Ei) obtained from 10 expo- 
nential fits of the form in equation [5] to a 3 x 3 matrix of T 
correlators as described in the text. The 4th column gives 
the T mass as determined from eq. [7| Fewer configurations 
were used for this than for the full calculation (and given in 
Table [TTJ) in several cases. For set 3, 202 configurations were 
used and for set 8, 470 configurations. The 5th column gives 
the result for the lattice spacing from setting the 2S — IS 
splitting equal to the experimental value of 0.5630 GeV |15| . 



Set 


aE 1 


a(E 2 - £i) 


aM 


a T /fm 


1 

2 


0.28775(8) 
0.28814(8) 


0.4244(33) 
0.4309(32) 


7.226(12) 
7.231(12) 


0.1488(12) 
0.1510(11) 


3 
4 


0.29330(3) 
0.29261(6) 


0.3439(8) 
0.3462(38) 


5.983(10) 
5.985(11) 


0.1205(3) 
0.1213(13) 


6 


0.26618(5) 


0.2381(37) 


4.281(12) 


0.0835(13) 


8 


0.24850(3) 


0.1679(14) 


3.050(18) 


0.0588(5) 



were estimated in [7] at 0.7% on the coarse and 0.6% on 
the fine lattices, so we include an overall systematic error 
of 0.7% in our error analysis here. 

One ingredient missing from our calculation and 
present in the experimental world is electromagnetism. 
This is then another possible source of systematic error. 
From a potential model calculation we estimate the shift 
in the 2S — 15 splitting to be less than 1 MeV from the 
Coulomb interaction between b and b (the electromag- 
netic self-interaction is included in the b quark mass). 
At less than 0.2%, this is negligible. 

To extract a physical value for the static-quark poten- 
tial param eter ri , we must combine the lattice spacings 
aj in Table " 



III 



. with the corresponding values of {r\/a)i in 
Table [T] and extrapolate to zero lattice spacing, correct- 
ing the sea-quark masses. We do this by fitting (ri/a)iaj 
from the ith ensemble to a formula for the effective ri 
corresponding to mx' — ray: 



rJ(a,8mT !l ,5m7 a ) = n 



(8) 



x 1 + c B 



25mf a + Smf a 



where r\ (the extrapolated value), c sea and cj are the pa- 
rameters tuned by the fit. Here the Sm se& are differences 
between the sea-quark masses used in the simulation and 
the correct masses for I = u/d and s quarks (see Ap- 
pendix [C]) . 

We have included twice as many terms as we need in 
the expansion in ajr\\ taking half as many terms gives 
essentially identical results. We are able to retain higher- 
order terms because we include Bayesian priors in our 
fit for each expansion coefficient used — that is, we in- 
clude an initial estimate for each parameter. Each prior 
functions as an additional piece of input data in the fit, 
thereby guaranteeing that we always have more fit data 
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TABLE IV: Major sources of uncertainty in physical n values 
obtained from simulation results for m-r' —mr, rriD B —mri /2, 
f Va , and from an analysis that fits all three types of data 
simultaneously. 





T 


D s 


f,, 

JVs 


combined 


a 2 extrapolation 


0.4% 


0.8% 


0.1% 


0.2% 


m s extrapolations 




0.0 


0.0 


0.0 


ri/a uncertainty 


0.2 


0.2 


0.1 


0.1 


initial uncertainty in n 


0.6 


0.8 






TY-K-r/g analysis 






0.8 


0.6 


statistical errors 


1.1 


0.4 


0.2 


0.3 


sea-quark mass tuning 


0.1 


0.2 


0.0 


0.1 


overall systematic error 


0.6 


1.1 




0.2 


Total 


1.4% 


1.6% 


0.9% 


0.7% 



than parameters, no matter how many parameters we 
choose to keep. In Bayesian fits like ours it is important 
to keep more parameters rather than fewer, not because 
they improve the fit but rather because they help us avoid 
underestimating our extrapolation errors. Here we used 
priors cj =0(1), cj ca — 0.0(1), both of which are broader 
(i.e., more conservative) than suggested by the empirical 
Bayes criterion |14j . Setting cj ca to zero has negligible 
impact, so this parameter is not really necessary. We 
also take a very broad prior for r x that encompasses all 
current estimates: r\ — 0.315(10) fm. It has little impact 
on our final errors. 

Fitting our data, we obtain a final value for the phys- 
ical T\ from the upsilon simulations of: 



ri = 0.3091(44) fm (from m T < - m T ). 



(9) 



The fit is excellent, with a % 2 per degree of freedom of 0.2. 
We show plots in Section III. The main sources of error in 
this result are listed in the T-column of Table IV most 
of the error is due to statistical errors from the Monte 
Carlo simulation. Improvement will require much higher 
statistics on the fine and superfine lattices. 



B. 



- m Vc /2 



A useful mass difference in the charm sector is that 
between the D s meson and one half of the mass of a 
low-lying cc state. We choose the r\ c rather than the 
J/ip because it is the easiest for us to calculate. This 
splitting has the experimental value 0.4784(7) GeV [TS] 
which changes to 0.672(2) GeV when the charm quark is 
replaced by a bottom quark. So the sensitivity to the 
heavy-quark mass, while stronger than for the heavy- 
onium splittings, is still rather mild. Both the D s and 
the rj c are ground state mesons and do not have the poor 
signal/noise issues that the T(2S') state had in the pre- 
vious subsection. Two quark masses are involved in the 
D s — ri c /2 splitting, however, and this makes the tuning 
rather complicated. As a result we have only done this on 
two ensembles - coarse set 4 (see Table ^ using 595 con- 
figurations with 2 time sources for propagators on each 
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FIG. 2: Results for m Vc vs m VB for different charm and 
strange quark masses on the coarse 01/05 ensemble (set 4). 
Points corresponding to different charm quark masses are 
given in different colors. Several different strange quark 
masses are given for each charm quark mass. The lattice 
spacing is determined from mo, — m I)c /2. Note that the ex- 
perimental point is shifted to allow for electromagnetic effects 
missing from our calculation, as described in the text. 



configuration and fine set 6 using 566 configurations with 
4 time sources each. 

Following the development of the HISQ action [17] it 
is now possible to handle charm quarks accurately with 
a relativistic action in lattice QCD. The success of this 
action is demonstrated by an accurate (7MeV) determi- 
nation, that agrees with experiment, of the masses of D 
and D s meson when the charm quark mass is fixed from 
the r] c pp. This has not been possible with any other 
discretisation of QCD for charm quarks. Here we are es- 
sentially inverting this calculation to use raj s — m„ c /2 to 
determine a value for the lattice spacing, simultaneously 
requiring that the mass of the r\ c and the r\ a (the fictitious 
pseudoscalar particle made of an s quark-antiquark pair, 
see subsection C) be correct. We use the HISQ action 
as described in [T7] except that we simplify the tuning of 
the coefficient of the Naik term (that corrects for a 2 er- 
rors) so that it is correct as a function of ma at tree level. 
In |17j it was shown that a nonperturbative tuning of this 
coefficient gave results very similar to the tree-level re- 
sult and so it is much simpler to take the tree-level result 
at each value of the quark mass. The difference between 
tree-level and nonperturbative tuning of the coefficient is 
a small discretisation error at relatively high order, and 
it will be taken care of in our continuum extrapolation. 

We proceed by calculating r] c , r] s and D s correlators for 
several different combinations of bare quark masses for 
charm and strange quarks. We use random wall sources 
as for the b quarks in the previous subsection, which im- 
proves the statistical error on the ground state masses 
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that we need here significantly. No smearing function is 
necessary. Each of the correlators is fit to an appropri- 
ate multi-exponential form including oscillating states for 
the D s pQ. 



fe=i 



a k (e 



-M k t 



-M k (T-t) 



), 



(10) 



for rj c and rj s and 



G mcson (t) = ^(e-^ + e-W-*)) (11) 

k=l 
fco=l 

for D s . As before, standard Bayesian fitting techniques 
are used |14j and results are taken from fits with n exp = 4, 
where we find the results, and their errors, are stable to 
changes in n exp . 

For each combination it is then possible to plot the 
value of m I)c against that of m Vs using the lattice spacing 



from m £> s 



J 2 and interpolate to experiment. A plot 



showing our results on the coarse 01/05 (set 4) ensemble 
is given in Figure [2] Interpolation to the matching point 
is relatively simple because the D s mass is linear in m s a 
and m c a, the r\ c mass is linear in m c a and the square of 
the rj s mass is linear in m s a for small changes in the quark 
masses. As can be seen from the Figure it is possible to 
pinpoint the matching point precisely and then confirm 
it with an additional calculation. Statistical errors below 
0.5% are achievable. 

The experimental values to be used here must take 
into account that the lattice calculation is missing the 
electromagnetism of the real world. This gives a signifi- 
cant shift to mo, — ti?7 c /2 because the D s is electrically 
charged and the r\ c is neutral, so their masses shift in 
opposite directions and the shifts add together. We es- 
timate the electromagnetic shift to the D s mass to be 
1.4MeV, and the shift to the ry c mass to be -2.6 MeV. In 
addition the r\ c in the real world can annihilate to gluons 
but not in our calculation and we estimate the shift from 
this effect to be -2.4MeV [T7]. This shifts m Ds - m Va /2 
from the real world to 0.4745 MeV for comparison to our 
calculation. We take the systematic error on this value 
to be one half of the shift we apply: that is, 2 MeV. This 
0.5% error gives a 1.5% systematic error on the values of 
the lattice spacing that we obtain, unfortunately domi- 
nating our statistical errors. The experimental mass for 
the rj c becomes 2.985 GeV after applying the shifts above 
and the r] s mass is taken as 0.686 GeV (see Appendix |A| 
and Eq. ([A6| ). 

Table |V| gives results on the coarse (set 4) and fine (set 
6) lattices from this method. The i] s correlators are a 
subset of those used in subsection C. The results for m, h 
are very slightly different in the two cases because of the 
different fitting strategy employed. 

Given the lattice spacings in Table |Vj we can again 
combine them with the corresponding values for r\/a 



TABLE V: Results for A = m_o s — m Vc /2 in lattice units from 
different charm and strange HISQ quark masses on ensembles 
4 and 6. The corresponding lattice spacing, at the tuned point 
where m Vc and m Vs agree with their physical values, are given 
in the bottom row of each section of the table. Errors shown 
are from statistics and extrapolation only. 

set 4 



m c a m Vc a 

0.72 1.98114(15) 

0.753 2.04293(10) 

0.753 2.04293(10) 



m s a m Vs a 

0.06 0.45787(23) 

0.06 0.45787(23) 

0.063 0.46937(24) 



aA 
0.3180(5) 
0.3214(5) 
0.3247(5) 



a/fm 



tuned 



0.3247(5) 0.1350(2) 



set 6 



m c a m Vc a m s a m 7]s a aA 

0.44 1.33816(7) 0.0358 0.30332(12) 0.21244(23) 

0.44 1.33816(7) 0.0382 0.31362(14) 0.21535(22) 

0.45 1.35934(7) 0.0358 0.30332(12) 0.21350(24) 

0.45 1.35934(7) 0.0382 0.31362(14) 0.21640(23) 



a/fm 



tuned 



0.2174(5) 0.0904(2) 



from Table [i] to obtain effective values rf s that we can 
extrapolate to zero lattice spacing. We do this using the 
same parameterization and priors for r^ s as we did for 
T\ in the previous section, except that here we allow for 
less dependence on the sea-quark masses since that is 
what our previous simulations have shown PQ. We take 
c sca = 0.00(1) as a prior. With only two data points, the 
fit is almost trivial, giving 



ri = 0.3157(53) fm (from mr>. 



= /2), (12) 



which agrees well with our estimate from the T but is 
less accurate. The main sources of error in this result are 
listed in the D s -column of Table IV the largest source 
of error is the overall systematic error 1 18] . The system- 
atic error could be improved slightly by using the J/tp 
instead of the r\ c to avoid the sizeable mass shift from 
annihilation to gluons and its uncertainty. In addition 
a more accurate understanding of electromagnetic mass 
shifts, with quantitative tests on the lattice would help 
(see, for example, [T§]). 



fv. 



The rj s is a fictitious pseudoscalar meson. It is like the 
pion and kaon, but with valence ss quarks. In the real 
world the valence ss state mixes with uu and dd, through 
valence quark- antiquark annihilation, to form the r\ and 
rf mesons. By omitting valence quark-antiquark anni- 
hilation from our simulation, we obtain the rj s instead. 
This meson is easily studied, in lattice QCD, using sim- 
ulations and, in the continuum, using partially-quenched 
chiral perturbation theory [5U]. In Appendix |A| we show 
how to determine its mass and decay constant from sim- 
ulation and experimental data for pions and kaons, using 
chiral perturbation theory. We are able to determine 
both parameters to within about 0.5%. 
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TABLE VI: Simulation results for the rj s mass m,, and decay 
constant f Vs for several lattice parameter sets (see Table [l| 
and s-quark masses am a . We also list the number of gauge 
field configurations and time sources per configuration used. 



Set 


am 3 


af Vs 


am Vs 


n c f g X nt 


1 


0.066 


0.1429(4) 


0.5250(6) 


631 x 2 




0.08 


0.1485(4) 


0.5782(6) 


631 x 2 


2 


0.066 


0.1436(4) 


0.5248(6) 


631 x 2 


3 


0.0537 


0.1144(2) 


0.4310(4) 


518 x 2 


4 


0.0546 


0.1160(3) 


0.4367(5) 


595 x 2 




0.05465 


0.1160(3) 


0.4369(5) 


595 x 2 




0.06 


0.1182(4) 


0.4580(5) 


595 x 2 


5 


0.0525 


0.1149(4) 


0.4259(6) 


460 x 2 




0.0556 


0.1161(4) 


0.4384(6) 


460 x 2 


6 


0.0358 


0.0806(2) 


0.3035(3) 


566 x 4 




0.0366 


0.0810(2) 


0.3069(3) 


566 x 4 




0.0382 


0.0817(2) 


0.3137(3) 


566 x 4 


7 


0.03635 


0.0811(2) 


0.3050(4) 


265 x 4 


8 


0.024 


0.0556(1) 


0.2120(2) 


218 x 4 


9 


0.0165 


0.0408(1) 


0.1548(1) 


200 x 2 




0.018 


0.0417(2) 


0.1621(2) 


101 x 1 



Given an accurate physical value, the rj s mass is the 
easiest quantity to use for tuning the s-quark mass in 
lattice simulations. It is significantly simpler to use than 
the K mass since m rh , unlike rriK, is only weakly de- 
pendent upon the u/d mass and therefore requires only 
minimal chiral extrapolation. This is because u/d quarks 
enter only in the sea for this meson. Another advantage 
of the rj s is that it is much less expensive to simulate than 
the K. 

Given a tuned s-quark mass, f„ s is much more use- 
ful for tuning the lattice spacing than either fx or f n . 
Again, this is because it is almost independent of the u/d 
mass (and because it is much less expensive to compute). 
We have computed both the decay constant and the mass 
for the rj s for a variety of s-quark masses for all of our 
lattice parameter sets 



The results are given in Table VI 
As discussed earlier, we used the HISQ formalism for 
the valence s quarks in our analysis, together with the 
MILC gluon configurations described in Table |TJ We an- 
alyzed the rj s created by the partially conserved axial- 
vector current in the HISQ formalism, so that the de- 
cay constant is automatically correctly normalized, with 
no need for further renormalization constants. We used 
random-wall sources when computing quark propagators, 
as described earlier for &-quark propagators, and used 
sources on several time slices for each configuration, to 
increase statistics, see Table |VT| 

We extracted masses and decay constants from the me- 
son correlators by fitting the middle 40% of the t range to 
a single exponential. This is less sophisticated than our 
approach to fitting correlators in previous subsections, 
but it simplifies the analysis of statistical correlations 
between different results coming from the same ensemble 
(with different s-quark masses). We get identical results 
if we use instead results from multi-exponential fits, ig- 



noring correlations. The x 2 P er degree of freedom of our 
fits was larger than one for some ensembles, possibly be- 
cause of lower statistics. To be conservative, we doubled 
the sta tistical errors everywhere (giving the results in Ta- 
ble VI I, resulting in excellent x 2 s. 

In analyzing our simulation results for (a/^J.; and 
(am Vs )i, we need to account for three systematic effects. 
First none of the simulations has precisely the correct s- 
quark mass m s ■ We did simulations at multiple values of 
m s so that we could interpolate. The lattice spacing can- 
cels out in the ratio (af Va )i/ (am Vs )i; we in effect vary m s 
until this ratio has the correct continuum value, obtained 
from our chiral analysis (see Appendix |A| . 

The second important systematic effect is that our sim- 
ulations have finite-lattice-spacing errors. We model de- 
pendence on the lattice spacing using a power series in 
(a/ri) 2 . A final, but much less important systematic 
is that the sea-quark masses are not quite right in our 
simulations. We did simulations using several different 
sea-quark masses so that we could correct for this depen- 
dence, which, as discussed above, we expect (and find) 
to be very small. Other systematic errors are negligible. 
In particular, finite-volume corrections for the ry s are no 
larger than 0.1% in our simulations. 

We account for these systematic effects by fitting our 
results (af Vs )i from the simulation using ensemble set i, 
with each s mass, to: 

{a/n^rl'f^iaux^), (13) 

where again values for (?"i/a)j come from Table [TJ This 
formula defines r^ 3 , which is the effective value of im- 
plied (for each ensemble set) by our data for the rj s decay 
constant and mass. We parameterize r^ s the same way 
we parameterized r^f and r®': 



' (a, Srri' 



f a ,(5m soa 



)=n 



(14) 



25r 



h 



where again r x is the physical value. Function f^(a, x Vb ) 
models the s-quark mass dependence of the decay con- 
stant where 



(15) 



is a measure of difference between the correct s mass 
and the s mass used in the simulation to produce (af Va )i 
and (am, ls )i. We parameterize as follows: 



fl 1 f(a,x Vs )=f ris +J2dkX k rh 



(16) 



k=l 



We allow the first two terms in the expansion to depend 
upon the lattice spacing by taking 

dk = d kQ + d k i(a/ri) 2 (17) 
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for j — 1,2; lattice-spacing dependence in the higher- 
order terms would have negligible effect (as do the higher- 
order terms themselves, as it turns out). 

Again we have included twice as many terms as we 
need in the expansions in a/r\ and x Va ; taking half as 
many terms in both cases gives essentially identical re- 
sults. Here we used priors cj s = 0(1) and c^; a = 0.0(1), 
as before, and dki = 0.0(5). Again all priors are some- 
what broader (that is, more conservative) than suggested 
by the empirical Bayes criterion |14j . 

The other parameter varied in the fit is the contin- 



uum/physical n in Eq. (14 1. We tried two different priors 
for this parameter. First we took the very broad prior, 
0.315(10) fm, we used for the other quantities. We also fit 
using the r% result from our chiral analysis of and fx 
in Appendix [A| (Eq. ( |A6| ). These two choices give results 
that differ by only a tenth of a standard deviation, which 
is negligible. We use the latter choice for our results be- 
low. We also take the values for f Vs and (/t; 3 /to?? s ) used 
in Eqs. ( 16 1 and ( 15 1 from our chiral analysis as described 



in Appendix |A[ 

Our final result for the continuum value for r\ in this 
analysis is: 



= 0.3148(28) (5) fm (from / n J, 



(18) 



where, as discussed in Appendix [A] the second error cor- 
responds to uncertainty about finite-volume corrections 
in the chiral analysis. The fit is excellent, with a x 2 per 
degree of freedom of 0.4. The main sources of error in 



this result are listed in the f Vg -column of Table IV the 
largest source of error is uncertainty in the physical val- 
ues of f ris and m ris from the tt-K-t] s chiral analysis. 



III. TWO RECIPES 

Two accurate recipes for setting the lattice spacing fol- 
low from the analysis in the previous section. The first 
requires that the static-quark potential be computed in 
the simulation, and a value for r\ja extracted from the 
results. This has been done accurately by the MILC col- 
laboration for their ensembles and we use their numbers. 
T\ja can then be converted to a value for the lattice spac- 
ing by dividing into the physical value of r\ . In the pre- 
vious section, we did separate determinations of T\ using 
simulation results for the upsilon and D s mass splittings, 
and for the rj s decay constant. For each we extracted ef- 
fective values of r\ for each lattice ensemble and param- 
eter set; and we extrapolated to the continuum to obtain 
physical values for r\. We have also done a joint analysis 
of all three sets of simulation results which is identical 
to what we did for each separately, but requiring that 
each fit use the same physical t\ — that is, we require all 
three to agree on the final value for v\ . This analysis also 
implicitly includes the r\ result from our chiral analysis 
of fn and fx since we use that value as the input prior 
for the combined analysis. When we do this we obtain 
the following final result, where again the second error 




0.00 0.05 



0.10 0.15 

(«Ai) 2 



0.20 



FIG. 3: Simulation results for the effective ri obtained from 
mo B — m Vc /2 (top), f Vs (middle), and ra^ — mi (bottom) 
are plotted versus (a/ri) 2 for various values of the sea-quark 
mass. The lines show the tuned fit functions from our simul- 
taneous fit to all three sets of simulation results. We used 
the fit functions to correct the simulation data points for the 
sea-quark masses; data points and lines are for 5m S q a = 0. 
The gray band is the continuum value obtained from the fit: 
n = 0.3133(23) fm. 



is due to uncertainties in finite- volume corrections to the 
chiral analysis (see Appendix [A| : 



r x = 0.3133(23) (3) fm (combined) 



(19) 



The fit is excellent with a \ 2 P cr degree of freedom of 0.4. 
Figure [3] shows r\ values from all three simulations plot- 
ted against the square of the lattice spacing. The sources 
of error in this combi ned analysis are summarized in the 
last column of Table IV The a 2 dependence in that fig- 



ure is all relative to a 2 dependence in the r\ /a values ob- 
tained from the static-quark potential. Thus the upsilon 
analysis has a 2 errors most similar to those in the static- 
quark potential's r±/a, while the D s analysis has errors 
least like those coming from the static-quark potential. 
There is no way to tell which of these quantities has the 
smallest absolute finite-a errors from just this simulation 
data; all that we can say is that they are consistent with 
each other in the continuum limit. 

The second recipe for determining the lattice spacing 
for a particular configuration set requires only the eval- 
uation of the mass and decay constant for the rj s (see 
Section II C I on those configurations; there is no need for 
the static-quark potential in this recipe. Lattice results 
for a/ I?a are fit to the formula 



af v f = af Vs (l + c\x v , + c 2 x 2 s ) , 
where, as before, 

tl at \ 2 

-1, 



(20) 



(21) 



and a, c 1; and c-i are fit parameters. Physical values 
for the mass and decay constant, m ria and f ris , are again 
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taken from Eqs. ( A6 1 in Appendix|A] Our simulations in- 
dicate that ci = 0.33(5) and C2 = 0.0(5) are good priors 
for these parameters; further terms in the x rjs expansion 
are unnecessary provided x Vs is small (it is less than 0.06 
for our data). The lattice spacing for the particular con- 
figuration set under study is then an output from the 
fit. 

In a typical simulation one guesses a value for the bare 
s-quark mass, am s in lattice units, to use in the quark 
action. Provided this is close enough to the correct value, 
a fit of the rj s results from this single mass is enough to 
generate an accurate lattice spacing. Doing simulations 
with two or more s-quark masses improves the result. 

The correct value for am s can also be estimated using 
a formula similar to Eq. ( 20 ) . A simpler procedure that 



gives almost identical results (to within 1/4%) for the 
correctly tuned s mass uses 



TABLE VII: Lattice spacings (in fm) and s-quark masses (in 
lattice units) determined using our n and f Vs recipes. Results 
are given for each configuration set from Table [T] We also list 
the number of am s values used in the r] s recipe. Note that 
the estimates converge as the lattice spacings vanish. 



Set 


#am s 


a\ ri 


ami uned \ ri 


a \vs 


am t s uned |„ s 


1 


2 


0.1456(11) 


0.0613(12) 


0.1583(13) 


0.0724(15) 


2 


1 


0.1465(11) 


0.0622(12) 


0.1595(14) 


0.0736(16) 


3 


1 


0.1184(9) 


0.0489(9) 


0.1247(10) 


0.0542(11) 


4 


3 


0.1197(9) 


0.0495(9) 


0.1264(11) 


0.0553(11) 


5 


2 


0.1185(9) 


0.0491(9) 


0.1263(11) 


0.0558(12) 


6 


3 


0.0847(6) 


0.0337(6) 


0.0878(7) 


0.0362(7) 


7 


1 


0.0844(6) 


0.0336(6) 


0.0884(7) 


0.0369(7) 


8 


1 


0.0592(4) 


0.0226(4) 


0.0601(5) 


0.0233(5) 


9 


2 


0.0440(3) 


0.0161(3) 


0.0443(4) 


0.0163(3) 



..tuned 



" 7/ 

ilat 



(22) 



where the lattice spacing is obtained from one of the two 
recipes above (or any other). Again in typical simula- 
tions, one guesses a value for am s and then uses r) s re- 
sults for this mass, together with this formula, to refine 
the initial guess. 

We compare lattice spacings determined using each of 
our two recipes in Table |VII| As expected, the lattice 
spacings are very different on the coarser lattices. This 
is because a 2 errors differ between the r\ and r\ s measure- 
ments. Also as expected (and required), the two recipes 
converge for smaller lattice spacings, as a 2 errors in both 
types of measurement become negligible. The errors in 
each case are comparable. We also include values for the 
correctly tuned s-quark mass (in the HISQ formalism) for 
each configuration set, and for each recipe for the lattice 
spacing. 

In neither of our recipes do we attempt to correct for 
sea-quark masses that are not correctly tuned. This is 
standard practice in lattice determinations of the lattice 
spacing. It pushes any sea-quark mass dependence from 
r i or f rjs (or whatever is used to determine the lattice 
spacing) into the other measurements of interest. This is 
a small effect for r\ and , and it is typically extrap- 
olated away together with the sea-quark effects intrinsic 
to the other measurements. 



IV. r 

ro/a is not determined directly by the MILC Collabo- 
ration. Instead they determine the coefficient of the 1/r 
term in the static potential in the region 0.2 — 0.7 fm. 
If this coefficient is B then: 



r„ = B + C ro 

n V b + c ri 



(23) 



where C ro = 1.65 and C ri = 1.0 [3J. This assumes 
that the same constant 1 jr coefficient would be obtained 
around r w r Q and r ~ and there will be a small 
systematic error, yet to be determined [5] for this as- 
sumption. B shows dependence on the lattice spacing 
and the sea quark masses as demonstrated in Figure 13 
of 5 . Extrapolating to the continuum and chiral limits 
gives B = —0.464(7), implying from equation 23 with 
the caveats above, that r^/ri = 1.488(5). Our value for 
r\ then gives tq = 0.4661(38) fm. This is in agreement 
with, but more accurate than, the previous MILC deter- 
mination of 0.462(12) fm which used ensembles at fewer 
values of the lattice spacing, but which includes a sys- 
tematic error of 0.004 from the variation of results with 
fit range in r. Our result also agrees with the direct de- 
termination from Aoki et al [5T] of r$ =0.48(1)(1) fm, 
which also includes the effect of u, d and s sea quarks 
and comes from an analysis with multiple values of the 
lattice spacing. 



V. CONCLUSIONS 

The accurate determination of the lattice spacing is 
of critical importance to obtaining accurate results from 
lattice QCD. Here we give two ways to do this with sub- 
1% errors for the first time. 

The first method makes use of the ~ 0.3% accurate 
values for ri/a calculated by the MILC collaboration on 
their ensembles (which could also be reproduced on other 
ensembles with similar statistics) coupled with the 0.8% 
accurate value for r\ given here : v\ — 0.3133(23) fm. 
Our result is 1.5er from our previous analysis [7] using 
only the T IS — IS splitting on fewer ensembles and 
combined with less accurate r\ja values. It is also la 
lower than that of the MILC collaboration using essen- 
tially the same results [3J. It is in agreement with, but 
slightly more accurate than a newer result from MILC [5 
of n = 0.3108(+^°) fm using fn data across a similar 
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range of lattice spacing values to our / I)s analysis but with 
ASQTAD valence quarks rather than HISQ quarks [5]. 

The second method is possibly simpler (in the absence 
of r\ja values) since it relies only on a standard me- 
son spectrum calculation that would automatically be 
included in many lattice analyses. The mass and the 
decay constant of the r\ s can be determined to better 
than 0.25% given similar statistics to those we have used 
here and provided a quark formalism is used in which 
the PCAC relation holds so that the decay constant has 
no renormalisation. Then the physical values for f Vs and 
m, lg that have been determined here can be used to find 
both the tuned value of the strange mass, by interpo- 
lation in f, lB /m Vs to 0.2647(18) and the lattice spacing, 
from taking f rjg =0.1815(10) GeV at the tuned point. 

The two methods are compared for the MILC ensem- 
bles in Table PVTTl 

To improve these methods so that errors below 0.5% 
are possible will require improvements in the chiral anal- 
ysis determining the r/ s parameters. These can be gauged 
from the error budgets in Tables |TV| and |rX) Key improve- 
ments that are certainly possible are statistical errors 
in the lattice results and accurate lattice data closer to 
the chiral and continuum limits. Improvements to other 
methods of determining the lattice spacing, such as that 
using the T spectrum and mr> s — m I)c /2 discussed here 
are important for cross-checks of systematic effects. 
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FIG. 4: The pseudoscalar decay constants plotted versus 



quark mass; ml / (2r 



m^) is approximately the ratio of 



the u/d to s quark masses: mi/m 3 . The fit data is from lat- 
tice simulations with three different lattice spacings; results 
decrease with decreasing lattice spacing. The data have been 
adjusted to correspond to points where the sea-quark masses 
correspond to the valence masses. The lines are from the 
tuned fit function for each of the three lattice spacings. The 
bottom line in each group is the extrapolation to a = 0. The 
gray bands indicate final values from the fit for the physical 
decay constants for all three mesons; the leftmost data points 
for f-n and /k are the current experimental values. 



TABLE VIII: Simulation results for pseudoscalar meson de- 
cay constants and masses (in lattice units) for several different 
lattice parameter sets (see Table]!]), u/d valence-quark masses 
mi, and s valence-quark masses m s . 



Set ami 



an is 



air i.k 



1 0.0132 0.066 0.1152(3) 0.2408(6) 0.1290(4) 0.4081(6) 

2 0.0264 0.066 0.1254(4) 0.3348(6) 0.1345(4) 0.4399(7) 

3 0.0067 0.0537 0.0889(3) 0.1567(4) 0.1020(3) 0.3242(5) 

4 0.01365 0.05465 0.0957(4) 0.2222(5) 0.1060(4) 0.3463(6) 

5 0.0278 0.0525 0.1041(3) 0.3113(6) 0.1095(4) 0.3727(6) 

6 0.00705 0.0366 0.0647(2) 0.1377(4) 0.0731(3) 0.2375(4) 

7 0.01635 0.03635 0.0710(2) 0.2050(4) 0.0759(2) 0.2594(4) 



APPENDIX A: / w , f K AND f Vs 

In [TJ , we analyzed simulation results for pion and kaon 
masses and decay constants obtained using the HISQ ac- 
tion for the valence quarks, with gluon configurations 
from MILC, produced using the ASQTAD action for the 
(rif = 3) light sea quarks. We described how to extrap- 
olate these results to the correct light-quark masses and 
to zero lattice spacing, obtaining decay constants that 
agree well with experiment. 

Here we reuse our earlier simulation results, which arc 



summarized in Table VIII (and in Table VI for the rj s ), 
to extract a value for the static-quark potential parame- 
ter r\ . More importantly, we also extract from this anal- 
ysis continuum values for the mass and decay constant 



of the r\ s meson. The masses and decay constants in 
Table |VIII| are obtained using the procedure described 
in Section II C| for analyzing r\ s correlators; we treat all 
mesons the same way. 

To extract a continuum value for r%, we fit the decay 
constant data for lattice ensemble i in Table VIII| (and 



Table VI for af Vs ) to 



(Al) 



where / ps is the formula from Appendix pj and (a,b) la- 
bels the valence quarks: (1,1) for pions, (7, s) for kaons 
and (s,s) for r] s s. The mass parameters x a , Xb ■ ■ ■ are 
computed from the simulation masses in Table |VIII| Pa- 
rameter r\ enters through the lattice spacing, which we 
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take to be 



(A2) 



where values for {r\/a)i are given in Table [!] We fit data 
from the pion, kaon and rj s simultaneously since all of the 
fitting parameters are universal. 

Our analysis here differs in three ways from our pre- 
vious paper pQ. First we are including the ry s in our si- 
multaneous analysis of the different pseudoscalar mesons; 
before we only included n and K mesons. Second we 
have re-expressed chiral perturbation theory in terms of 
pion and kaon masses rather than quark masses. This 
simplifies the analysis and also gives more reliable esti- 
mates for infrared quantities like chiral logarithms. We 
take the pion and kaon masses corresponding to the sea 
quark masses from [3] for ensemble sets (3,4,6,7). Re- 
sults for the other ensembles are not published so we 
generate approximate meson masses to go with the sea 
quarks by mu ltiplyi ng the meson masses for the valence 
quarks (Table Villi by (m sca /m val ) 1/2 (after converting 
HISQ quark masses into ASQTAD quark masses using 
Eq. (CI)). Replacing quark masses with meson masses 
in the chiral formulas gives results that agree well with 
our previous results. 

The third difference from our earlier analysis is that 
here we require the fitting function to also fit experimen- 
tal results for and Jk at zero lattice spacing. We do 
this by treating the physical results as additional data 
to be fit, together with the simulation results, to a sin- 
gle parameterization. In our previous study we fit only 
simulation results, showing that these agreed with exper- 
imental data. Here our goal is different, as we seek an 
accurate value for r\. That value is the one that allows 
the same chiral formulas to fit both our lattice results 
and the experimental results; r\ is determined, in effect, 
from the experimental values for f T and Jk ■ 

Our simulations omit both electromagnetic corrections 
and isospin-breaking effects. Following [21], we remove 
leading-order errors of both sorts by using 



to 2 - = 1 

m K 2 



m 2 K o 



L K+ 



(A3) 

(1 + A b )(to 2 +-to 2 „)) (A4) 



for the physical masses of the pion and kaon. pa- 
ramctcriscs the violation of Dashen's Theorem which, in 
the chiral limit, states that the K + and n + have equal 
electromagnetic corrections, while the tt° and K° have 
none. We take = 1(1). Electromagnetic corrections 
are also removed from the standard definition of the de- 
cay constants, whose values we take to be [15] : 

U = 0.1304(5) GeV f K = 0.1555(9) GeV. (A5) 

The fitting parameters that are varied in the fit include 
all of the parameters that define / ps (see Appendix [b]) , 
as well as r\ itself. As discussed in Appendix [5] all pa- 
rameters have priors in our fits. For r\ we take a very 
broad prior, r\ — 0.315(10) fm, that easily encompasses 



TABLE IX: Extrapolation and other errors in our results 
from the chiral analysis of 7r, K, and r] s masses and decay 
constants. Finite-volume errors are dealt with separately (see 
text). 





n 




m Vs 


fvs /mys 


a 2 extrapolation 


0.6% 


0.2% 


0.3% 


0.3% 


m q extrapolations 


0.7 


0.2 


0.2 


0.2 


ri/a uncertainty 


0.2 


0.1 


0.1 


0.1 


initial uncertainty in n 


0.6 


0.1 


0.0 


0.1 


experimental errors in n, K 


0.2 


0.2 


0.2 


0.3 


statistical errors 


0.7 


0.3 


0.4 


0.5 


Total 


1.4% 


0.5% 


0.5% 


0.7% 



all current estimates; it has little impact on the final re- 
sults. 

The results of our fit are show in Figure [4] The fit is 
excellent, with a x 2 P er degree of freedom of 0.4. Our 
main results are physical (i.e., continuum) values for ri 
and for the decay constant and mass of the rj s : 

ri = 0.3190(45) (20) fm, (A6) 
f Va = 0.1815(10)(2) GeV, 
m„ s = 0.6858(38)(12)GeV 
f v Jm Vs = 0.2647(18) (1) 

By "physical" we mean extrapolated to zero lattice spac- 
ing and the correct, physical values for the quark masses. 
We quote two errors here. The first is the fitting error, 
representing uncertainties from simulation statistics, and 
from the chiral and lattice-spacing extrapolation. A de- 



tailed breakdown of these errors is given in Table IX 
The second error is equal to the size of the finite-volume 
correction. As discussed in Appendix B, finite- volume 
corrections are somewhat ambiguous for staggered-quark 
formalisms like HISQ. We choose to include finite- volume 
corrections, but, to be conservative, take half the size of 
the correction as an uncertainty. 

We need the physical ij s results for our analysis in Sec- 
tion II C of the T] s decay constant L . These fit results 
have statistical correlations with each other, as well as 
with the output value of ri, the values of r\/a used in 
the fit, and the simulation results for af rjs (Table VI I. 



We used the fit here to compute means and a covariancc 
matrix for all of these quantities, and this is used as input 



data in the f rjs analysis of Section II C 

Note that the values for the ?y s mass and decay constant 
agree to better than a percent with the leading-order 
expectations from chiral perturbation theory: (2to^ — 
to 2 ) 1 / 2 and 2Jk — fn, respectively. Our analysis above, 
however, goes far beyond leading order (see Appendix[B]). 
Our r] s results are also quite independent of the input 
prior for r\\ taking 0.3133(23) fm as the prior, for ex- 
ample, causes shifts that are smaller than a quarter of 
a standard deviation. The rj s parameters are most sen- 
sitive to the physical parameters for the pion and kaon. 
They can easily be corrected should there be small shifts 
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in the values derived from experiment for /„• and /#-. The 
changes in the rj s parameters would be: 



A/„, - 0.6A/ K + 0.2AU 
Am,, = 0.2A/ K - 0.6 A f n , 



(A7) 
(A8) 



where Af n and Afx are changes in the pion and kaon 
decay constants from the values used here. 



APPENDIX B: AUGMENTED CHIRAL 
FORMULAS 

We model light-quark pseudoscalar masses and decay 
constants using partially-quenched chiral perturbation 
theory, augmented with corrections for the finite lattice 
spacing. For simplicity we re-express chiral perturbation 
theory in terms of pion and kaon masses, in place of the 
quark masses, using 



x, 



</2 



0.007 



ml/2 



0.17 



A x =47r/ ff /V2 « 1.2 GeV. 



A 2 

x 



as expansion parameters, where 



(Bl) 
(B2) 

(B3) 



We use the formulas through next-to-leading order 
from 20j, together with higher-order corrections in xi 
and x s and finite-a corrections. For example, we model 
the mass and lattice spacing dependence of the decay 
constants using 



,a) = f^ + Sfx + Sfi 



lat 



(B4) 



where J NLO is the chiral formula through next-to-leading 
order, 5f x is the continuum correction due to higher- 
order mass corrections, (5/i a t is the correction due to the 
finite lattice spacing, and (a, b) labels the valence quarks: 
(1,1) for pions, (7,s) for kaons, and (s,s) for rj s s. 

Our simulation results are not sufficiently accurate to 
resolve the difference between high-order polynomials in 
xi and x s and high-order logarithms, so we keep just the 
polynomials: 



(B5) 



8 fx = fo (ci(x a + x b ) 2 + c 2 (x a - x b ) 2 

+ c 3 (x a + x b )(2x s r+xr) 

+ c 4 (2xT a + xf a ) 2 + c 5 (2(xT a ) 2 + (xf a ) 2 ) 
+c 6 (x a + x b ) 3 + c 7 (x a + x b )(x a - x b ) 2 ) , 

where /o is the bare decay constant in chiral perturbation 
theory and the c« are expected to be 0(1), except for sea- 
quark terms where the coefficients should be 3-5 times 
smaller. Still higher-order terms are smaller than 0.1% 
and so negligible, as are the last few terms in practice. 
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FIG. 5: Fits to two different sets of fake data for pion and 
kaon decay constants with very different a 2 behavior from 
each other and from the real simulation data (Figure [4||. The 
"experimental" points indicated in each case correspond to 
the exact results, extracted from the formulas used to generate 
the fake data. 



Following [1J, we model the a 2 dependence using a mix- 
ture of terms that depend upon a 2 and xi : 



Sfut = fo (rfi(aA Q cD) 2 a s + d 2 (aA QC D) 2 a 3 
+ rf3(aA QC D) 2 log(a; / )a^ 
+d 4 (aA QCD ) 4 + d 5 (aA QCD ) £ 



(B6) 



where we set a s = ay(2/a), Aq CD = (Am^ — ml)/ '3 as 
in [5D]. (Aqcd is also the ultraviolet scale in the chi- 
ral logarithms.). We allow the coefficients to have mass 
dependence 

di = da + d i2 (x a + x b ) (B7) 

+d l3 (2x s r+x B r)+d l4 (x 2 a + x 2 ). 

Again the are expected to be O(l), except for terms 
involving sea-quark terms which should be 3-5 times 
smaller. The highest-order terms in these expansions are 
already negligible, making further terms irrelevant. We 
include the log(x/) term in Eq. (B6l to allow for non- 
analytic behavior at small xi, although in practice it is 
negligible in our fits. 

We included priors in our fitting analysis for each of the 
parameters in J NLO and for all the c.;S and dijS. These are 
initial estimates for each parameter that function as extra 



14 



"data" and allow us to account (in our error estimates) 
for the uncertainties in these parameters, even when they 
are largely unconstrained by our simulation data. The 
parameters in J NLO are well determined by our data; we 
use very broad priors for these, which have no impact on 
the final errors. We use a prior of 0(1) for each of the CjS 
and dijS, except for terms involving sea-quark masses in 
which case we use 0.0(3). 

As reported in [T], we have tested these fitting formu- 
las extensively by using formulas from partially-quenched 
staggered chiral perturbation theory, with randomly se- 
lected coefficients and randomly generated higher-order 
corrections in the masses and a 2 , to generate fake data 
sets for the same masses and lattice spacings used in our 
analysis here. We added statistical noise to the fake data 
that was comparable in magnitude to that in our real 
simulation data, with similar correlations. We then fit 
the fake data using the formulas above, together with 
the Empirical Bayes method [H] to set a prior for the 
expansion parameters (q and dij). In each case we could 
compare extrapolated results from our analysis of the 
fake data with the exact results, since we knew the un- 
derlying formula used to generate the fake data. We ran 
tests for several hundred cases. As expected, we found 
that 70% of the time the extrapolated results were within 
one standard deviation of the exact results. Two exam- 
ples, shown in Fig. [5] illustrate how effective our formulas 
are in handling a 2 dependence that is much larger and 
much more complex than we see in our actual simulation 
results (Figure [3J, 

The logarithms in the NLO chiral formulas reflect in- 
frared sensitivity. These terms are sensitive to the finite 
volume of our lattice at the level of 0.1-1% for the de- 
cay constants (less for masses). We add finite- volume 
corrections to the logarithms which we obtain by recom- 
puting the one-loop chiral corrections that lead to log- 
arithms using finite-volume sums instead of integrals in 
momentum-space, and subtracting them from the infi- 
nite volume results. These corrections are quite sensi- 
tive to the meson mass, which raises an issue since in 
staggered-quark formalisms like HISQ each pseudoscalar 
meson comes in several different "tastes" , all of them 
heavier than the Goldstone meson whose mass we use in 
our formulas. Taste splittings are a 2 corrections, which 
vanish in the continuum limit, and most of the effects 



we have tested, as discussed in the previous paragraph). 
The finite-volume corrections, however, are particularly 
sensitive to meson masses, so we use an "effective" pseu- 
doscalar mass when we calculate them: 



(m e M) 2 =(mf b r+g m (a/n) 



(B8) 



where m^ b is the Goldstone meson's mass. We expect 
g m « 0.2 GeV 2 . We allow g m to float in our fits, treating 
it as a fit parameter. We use 0.2(6) as our prior. Our fit 
favors a nonzero value for g m , giving g m = 0.2(3) which 
is consistent with expectations. 

APPENDIX C: SEA-QUARK MASSES 

We include terms in each of our fitting functions that 
correct for the discrepancies 5m s q c& between the bare sea- 
quark masses used in the simulation and the physically 
correct bare quark masses (that is, the ones that give cor- 
rect masses for the it, K, and r] s ). Our estimates for the 
correct s-quark masses (in lattice units) for each ensem- 
ble are given in Table VII| the u/d mass is 27.8(3) times 
smaller p]. These masses, however, are for HISQ quarks, 
while the sea quarks were all analyzed using the ASQ- 
TAD formalism. Quark masses in the two formalisms can 
be related to each other, ensemble by ensemble, by com- 
paring 7r and rj s masses for mesons whose valence quarks 
are either HISQ or ASQTAD quarks. HISQ masses and 
ASQTAD masses are equivalent when they give the same 
7r and rj s masses. The ratio of a HISQ mass to the corre- 
sponding ASQTAD mass determined in this way should 
be almost independent of the valence-quark mass, but 
will depend somewhat on the lattice spacing and weakly 
on the sea-quark masses. We have compared ASQTAD 
data from [3] for ensembles 3,4,6,7 with our results in 
Table |VHT| to obtain the following simplified parameter- 
ization for the ratio of HISQ to ASQTAD quark masses: 



am 



hisq 



flm asq 



= 1.158 



1 + 0.44 (a/ri) 2 



1 + 0.009 (am^/aml^) 



(CI) 



where 
and 



-.tuned 



asq 



is the tuned HISQ mass given in Table 



VII 



is the sum of the three sea-quark ASQTAD 



of these we model with our corrections Eq. (B6l (which 



masses for that ensemble, 
few percent. 



This formula is accurate to a 
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